00 

o 

O 
(N 

> 

O 

0^ 






I 

o 
o 



> 

m 

oo 

o 

oo 
O 



% 



Numerical evidence for unstable magnons at high fields in the Heisenber^ 
antiferromagnet on the square lattice 

Olav F. Syljuasen^ 

^Department of Physics, University of Oslo, P. O. Box 1048 Blindern, N-0316 Oslo, Norway 

(Dated: November 19, 2008) 

We find evidence for decaying magnons at strong magnetic field in the square lattice spin- 1/2 
Heisenberg antiferromagnet. The results are obtained using Quantum Monte Carlo simulations 
combined with a Bayesian inference technique to obtain dynamics and are consistent with predictions 
from spin wave theory. 
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The square lattice spin- 1/2 Heisenberg antiferromag- 
net (2DHAF) is the archetype quantum antiferromagnet 
and describes the magnetism of the mother compounds of 
the high-Tc cuprate materials. While the 2DHAF itself 
is rather well understood, less is known when additional 
interactions are added. In this Rapid Communication we 
report on the nature of the excitations of the 2DHAF in 
an external magnetic field. Such a perturbation is highly 
relevant from an experimental point of view and has very 
interesting consequences. 

Adding a magnetic field introduces interactions among 
the magnon excitations of the pure 2DHAF. According 
to spin wave theory these interactions will cause the 
magnons in certain regions of the Brillouin zone to be 
unstable at very high fields py]. One should however note 
that spin wave theory, which is an expansion in the pa- 
rameter 1/ S, is not at all guaranteed to work for S" = 1/2. 
Furthermore the spin wave prediction involves certain as- 
sumptions as it is a self-consistent calculation which ne- 
glects vertex corrections and renormalizes the spin wave 
spectrum considerably. In addition there are also other 
theoretical predictions of what happens to the 2DHAF 
in a magnetic field: A scenario where the spin-1 excita- 
tion is viewed as a composite excitation of two dimen- 
sional spinous predicts additional bands of low-energy 
excitations at moderate values of the field '4 [Sj. Thus 
to settle the issue an independent calculation of the ex- 
citation spectrum of the 2DIIAF in a magnetic field is 
needed. It is rather surprising that no such calculation 
has been performed until now, given the prominence and 
simplicity of the model. In this Rapid Communication wc 
use Quantum Monte Carlo (QMC) simulations combined 
with a Bayesian inference technique to give a detailed 
and unbiased description of the excitation spectrum of 
the 2DIIAF in a magnetic field. We find regions with 
broad spectral features at strong magnetic fields indicat- 
ing decaying magnons consistent with spin wave theory, 
and we also show details of the spectrum in the decaying 
regions. 

In zero magnetic field the ground state of the 2DHAF 
is Neel ordered with a renormalized moment. Spin wave 
theory explains very well the renormalized moment as 



well as the dispersion relation of the spin-1 excitations, 
except for an anomaly occurring at (7r,0) (we use units 
where the lattice spacing a = 1) where the dispersion 
softens and the spectrum broadens anomalously[J,l5|, ISJ. 
This anomaly at (tt, 0) is on the contrary natural in a pic- 
ture of the spin-1 excitation as a composite excitation of 
two spin 1/2 spinous that are excitations about a mean- 
field pi-flux state [2|- While this bound state has almost 
the same dispersion as predicted by spin wave theory it 
differs around (tt, 0) because the two spinous making up 
the bound state each has low energy close to the nodal 
pockets at (±7r/2, 7r/2). Furthermore this spinon picture 
can also explain the anomalous broadening of the contin- 
uum part of the spectrum at (tt, 0)[7|. Thus the spinon 
picture works well for explaining the (tt, 0) anomaly quali- 
tatively, however being essentially a projected mean field 
calculation it is not of the level of precision needed in 
order to compare quantitative details with exact results 
and conventional spin wave theory. 

In a magnetic field the distinction between the spinon 
and magnon picture is more dramatic. By assuming that 
the spinon bands move rigidly when changing the mag- 
netic field, as happens in the ID Heisenberg antiferro- 
magnet, the dynamic structure factor was calculated in 
Ref. [3( . The main result of this is the prediction of a new 
band of excitations forming an incommensurate ring at 
low energy close to (tt, tt) at intermediate fields. In con- 
trast, spin wave theory predicts a smooth change of the 
magnon dispersion with field up to a high value at which 
the magnons in some regions of the Brillouin zone become 
unstable and decay [l|. 

The Hamiltonian of the 2DHAF in a magnetic field is 



H = jY^S,-S,-BY^St 



(1) 
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where Si are spin-1/2 operators, (i, j) indicates nearest- 
neighbor lattice sites on a square lattice, J is the ex- 
change constant and B is the strength of the magnetic 
field, here pointing along the spin-z axis. 

This Hamiltonian can be simulated efficiently using 
QMC methods. Specifically we use the stochastic se- 
ries expansion technique with directed- loop updates [8|. 



In order to measure the excitations of the system we will 
focus on calculating the dynamic structure factor which 
is the primary quantity measured in neutron scattering 



(2) 
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S"iq,co)= / dte^-'{S^it)S%^{0)) 



where 5*2 = (1/A^) J2r^? and N ^ L'^ is the total num- 
ber of lattice sites. We use L = 32 when not stated 
otherwise. The spin component direction is denoted by 
a and can be either x,y or z. 

The dynamic structure factor is not directly accessi- 
ble in QMC because the method is formulated in imag- 
inary time. In contrast, what is accessible in the QMC 
simulations is the imaginary-time correlation function 
D"{q,T) = (5"(g,r)5'"(-(f,0)) which is related to the 
dynamic structure factor as 



n'i,.r,-l ^ 
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S"iq,Lo). (3) 



A direct inversion of Eq. ([3]) is difficult as the kernel has 
many near zero eigenvalues which causes the errors of 
D to be amplified enormously. Instead we use a statis- 
tical inversion method known as the Average Spectrum 
Method( ASM) 19|] where one takes as the resulting struc- 
ture factor the weighted average of all possible structure 
factors 5',"^ {q, uj) where the weight factor depends on how 
well the associated imaginary-time correlation function 
i??N matches the imaginary-time data D". In practice 
this weighted average is computed using another Monte 
Carlo simulation, and details are available in Ref . [10[ . 

Before we present the results for the dynamic struc- 
ture factor, we show in Fig. [T] the instantaneous spin- 
correlation function D"{q,0) for longitudinal (along the 
magnetic field) and transverse spin components, respec- 
tively. Using Eq. ^ these results translate into a mea- 
sure of the total strength of excitations at a given q- 
point. For the longitudinal spin component the effect 
of the magnetic field is to weaken the antiferromagnetic 
signal at (tt, tt) and enhance the ferromagnetic signal at 
(0, 0) reflecting uniform canting of spins along the mag- 
netic field. In contrast the antiferromagnetic signal at 
(tt, tt) remains even at high magnetic fields for the trans- 
verse spin components. For high fields note that D°'{q, 0) 
becomes nearly independent of q except close to q = 
and q= (tt, tt) for the longitudinal and transverse polar- 
ization respectively. 

In a weak magnetic field we expect that the dynamic 
structure factor is essentially the same as in zero field. 
Fig. [5] shows a color-scale plot of S^{q, cu) along the same 
scan in g-space as in Fig. [T] for B/J ~ 0.05 at T = J/40. 
The color-plot is made by discretizing the frequency in- 
terval uj/J = [0, 8] into 1000 bins and assigning a color to 
the value of S{q, uj) within each interval. If the value of 
S{q^Lo) corresponds to a value higher than the maximum 
color value, the color is set to the maximum color value. 
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FIG. 1: Instantaneous correlation function D°'{q,0) vs. q 
along the path in the Brillouin zone indicated in the upper left 
corner of the top panel. The different curves are for different 
values of the magnetic field B specified in the legends. The 
upper panel is for longitudinal spin polarizations (a — z) 
while the lower is for transverse polarizations {a = x). The 
error bars are smaller than the symbol size. 



In accordance with the zero- field result [5[ we find that 
the magnon energy is lower at (tt, 0) than at (7r/2,7r/2). 
More precisely, the location of the peaks in S{q, uj) at the 
two values of momenta are i-L;poak.(7r,o) = 2.21 ± 0.02 and 
^poak,(7r/2,7r/2) — 2.42±0.02, Corresponding to a 9% diflFer- 
ence. While the dispersing magnon peak is the dominant 
feature of Fig. [5] it is also possible to get a glimpse of 
high energy features. Noteworthy is the broadness of the 
magnon peak at (tt, 0) and the V-shaped feature above 
(tt, tt). 

In a weak magnetic field the spins will orient them- 
selves in a staggered configuration almost perpendicular 
to the field. As the field increases the spins will cant 
more and more along the field, and finally at B = AJ 
will align completely with the external magnetic field. 
From a symmetry perspective a magnetic field will break 
the SU(2) spin symmetry of the 2DIIAF. This breaking 
of a global symmetry is reflected in the spectrum by the 
fact that one of the two goldstone modes gets a mass. 
In the longitudinal channel this manifests itself at (tt, tt) 
where the zero-energy mode acquires a gap equal to the 
value of the magnetic fleld. In the transverse channel it 
is the q = mode that becomes gapped. Fig. [3] shows 
the longitudinal dynamic structure factor for B/J = 1. 



The dispersion seen in Fig. [3] agrees quantitatively with 
spin wave theory. By parameterizing the Hamihonian 
Eq. (IT|) in a canted coordinate system, characterized by 
the canting angle 9, and expressing the spin operators 
by boson operators according to the Holstein-PrimakofF 
transformation one gets terms with all orders of bosonic 
creation and annihilation operators, including linear and 
cubic terms. Minimizing the energy of the zeroth order 
term gives a condition relating the canting angle to the 
magnetic field: sinO = B/AJ . With this value of 9 the 
linear term vanishes identically, and the quadratic term 
can be diagonalized, yielding 



ojI = 2Jy/(l-7,,)(l + cos207£) 



(4) 



where 7^ = (cosfc^; + cosky)/2. There are corrections to 
this dispersion coming from higher order terms. The 3rd 
and 4th order terms in the expansion of the Hamiltonian 
are 



(«j> 



(5) 



H^'^' — — y^ < —2 cos29 TiiTij + sin'^ 9 a\ {ui + Uj) Uj 

- cos^ 9 a] ( flifli + a]a] j )aj + (i ^-> j) ^ (6) 

Following Ref. [l| we find corrections to the dispersion 
Eq. (HI) by calculating the self-energy S by contracting 
two cubic terms and retaining only the frequency inde- 
pendent parts of the quartic term. The perturbative cor- 
rection to the magnon dispersion is then 



^g=.;2 + E(a;"fc). 



(7) 



In zero field this reduces to a purely multiplicative renor- 
malization cjr = Zc^t^ where Zc = 1.1765 llj and is seen 
as a white curve in Fig. [2] while the white curve in Fig. [3] 
indicates ujj: according to Eq. ([7]) for _B = J. 

Fig. [3] reveals also that the broad peak seen at (tt, 0) 
is also present at B/J = 1. There are no appreciable 
continuum above (tt, tt), but high energy features are seen 
going away from (7r,7r). 

When the magnetic field is increased further our re- 
sults show that the dynamic structure factor is still dom- 
inated by a large dispersing magnon peak with a gap at 
(tt, tt) that continues to increase. This holds until the 
field reaches B « 3J, where the magnetization per site 
is M — 0.318, i.e. about 63% of the saturation value, at 
which the peaks at q-space points in the region around 
{tt/2, 7r/2) broadens and become less well-defined. Fig.lH 
which is the main result of this Rapid Communication, 
shows the dynamic structure factor at a field B/J = 3.5 
(M = 0.392) where these effects reach a maximum. 
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FIG. 2: Color-scale plot of S'^{Q,uj) where q takes values 
along a path in the Brillouin zone specified in the inset. The 
thin white curve indicates the renormalized linear spin wave 
dispersion, ZctOs^q)- 
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FIG. 3: Color-scale plot of S^ {q, uj) along the g*-path indicated 
in[T] The thin white curve is the perturbative result for the 
spin wave dispersion using the self-energy. 



For magnetic fields above iJ « 3J the self-energy S 
acquires an imaginary part due to decay of magnons thru 
the three-boson vertex in Eq. ([5]). This implies that the 
perturbation expression Eq. ([7]) will not work well for 
finding the dispersion and a more sophisticated analysis 
is needed [l|. An estimate for the imaginary-part of the 
self-energy, or equivalently the magnon linewidth, can 
however simply be obtained using Fermi's Golden rule 



where we use the bare spin-wave dispersion. The result 
of this is plotted in Fig. [1] and agrees qualitatively with 
the location of the decaying regions. This result depends 
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FIG. 4: Color-scale plot of S^ {q, lj) along the q*-patli indicated 
in Fig. [T]for B/J = 3.5. The thin black curve is the magnon 
linewidth F calculated using the golden rule Eq. ((S)). 
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FIG. 5: (Color online) Lineshapes of the dynamic structure 
factor at {n/2, tt/2) for three different values of magnetic field 
as indicated in the legend. Here the system size is L = 16 and 
the inverse temperature is /3J = 200. 



mostly on the available phase space, i.e. the spin-wave 
dispersion, as the interaction vertex is a slowly varying 
function of q. 

Note that the decay of a spin-1 excitation into two spin- 
1 excitations is not forbidden by the symmetry of the 
Hamiltonian. Not even at i? = where the Hamiltonian 
has a full SU(2) symmetry. However, at _B = there 
are no terms in the Hamiltonian with an odd number of 
boson terms, thus a single magnon can only decay into 
an odd number of magnons. In contrast cubic terms are 
present in a magnetic field. 

Fig. |4] reveals a double-peaked structure in the de- 
caying regions. Looking more closely at this we show 
in Fig. [5] lineshapes at {tt/2, 7r/2) for different values of 
the magnetic field. For the B/J = 3.5 curve one can 
clearly identify two peaks, one at lu ^ 1.28J and one 



at w sa 2.Q6J. However these peaks are overlapping. 
This could be either because the ASM cannot resolve the 
peaks completely or that the spectrum is really a con- 
tinuum where the two peaks indicate enhanced density 
of states near the edges. We favor the latter interpre- 
tation as the temperature is much less than the spacing 
between the peaks and we have collected very accurate 
QMC data (std.dev. ~ 1 x 10"*^). This should in our 
experience be enough to resolve the peaks if they were 
distinct [lO| . In view of the continuum interpretation the 
effect of increasing the magnetic field in this region is then 
to transfer weights from the lower to the upper edge of 
the continuum. 

In no circumstances do we see ring-like regions of low- 
energy excitations as predicted by the spinon bound state 
scenario. Thus our numerical results do not favor this 
interpretation. 

The magnon decay reported here is seen at very high 
magnetic fields. Experimentally this decay is therefore 
most relevant in materials with a low value of the ex- 
change constant [12l|. and perhaps also for systems with 
cold fermions in an optical lattice [13|. 

The simulations were carried out using the Titan clus- 
ter at the University of Oslo and the Nordugrid data 
production facility through the Nordugrid ARC middle- 
ware. 
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